!Entièrement validé, sauf:
!-routine de compression \lambda < 1; manque l'implémentation de la queue.
module fonctions
implicit none
integer, parameter :: points = 300 !nombre de points d'évaluation des fonctions
integer, parameter :: itmax = 80 !nombre maximal d'itérations
double precision, parameter :: infty = 10.d0 !valeur numérique pour l'infini
double precision, parameter :: lambda = 0.1d0 !vitesse de convergence (=1.0) versus stabilité (=0.01...) (?)
!Recollement OBLIGATOIRE si on veut résoudre pour c>1.
!L'équation pour petits c est utilisée jusque pour c=1, bande de recollement comprise.
double precision, parameter :: centrage = 0.9d0 !centrage du recollement < 1, largeur recoll: [cent - (1-cent); 1]
double precision, parameter :: pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164 
contains

!------------------------------------------------------------------------
!Méthode de simpson pour calculer une intégrale connaissant uniquement
!quelques points discrets. Explication des variables:
!  y       : tableau de dimension (n+1) contenant les valeurs y(xi)
!  h       : largeur du pas d'intégration
!  n       : nombre de pas d'intégration
!  somme   : valeur de l'intégrale sur la base des points fournis
subroutine simptab(y,h,n,somme)
implicit none
double precision, intent(in), dimension(0:points) :: y
double precision, intent(in) :: h
integer, intent(in) :: n
double precision, intent(out) :: somme
integer :: i
double precision :: wgt
somme = y(0)
somme = somme + y(n)
if (n==0) somme = 0.d0
if (n==1) somme = h*somme/2.d0
wgt = 2.d0
if (n>1) then
	do i=1,n-1
	  if (wgt == 2.d0) then
			wgt = 4.d0
	  else
			wgt = 2.d0
	  end if
		somme = somme + wgt*y(i)
	end do
	somme = h*somme/3.d0
end if
return
end subroutine simptab

!------------------------------------------------------------------------
!Méthode du trapèze pour calculer une intégrale connaissant uniquement
!quelques points discrets. Explication des variables:
!  y       : tableau de dimension (n+1) contenant les valeurs y(xi)
!  h       : largeur du pas d'intégration
!  n       : nombre de pas d'intégration
!  somme   : valeur de l'intégrale sur la base des points fournis
subroutine traptab(y,h,n,somme)
implicit none
double precision, intent(in), dimension(0:points) :: y
double precision, intent(in) :: h
integer, intent(in) :: n
double precision, intent(out) :: somme
integer :: i
somme = y(0)
somme = somme + y(n)
if (n==0) somme = 0.d0
if (n > 1) then
	do i=1,n-1
	  somme = somme + 2.d0*y(i)
	end do
end if
somme = h*somme/2.d0
return
end subroutine traptab

!------------------------------------------------------------------------
!Interpolation linéaire
!  a,b     : bornes d'interpolation, a < b.
!  fa,fb   : valeur de la fonction en ces points.
!  x       : point d'évaluation.
function IL(a,b,fa,fb,x)
double precision, intent(in) :: a,b,fa,fb,x
double precision :: IL
IL = (fb-fa)/(b-a)*(x-a) + fa
return
end function IL

!------------------------------------------------------------------------
!Normalisation de la solution, et de son second moment
!  y       : tableau de dimension (n+1) contenant les valeurs y(xi)
!  h       : largeur du pas d'intégration
!  l       : valeur de sortie du facteur de dilatation. A titre d'information.
subroutine normalise(y,h,l)
implicit none
double precision, intent(inout), dimension(0:points) :: y
double precision, intent(in) :: h
double precision, intent(inout) :: l !paramètre de dilatation lambda
double precision, dimension(0:points) :: tmp1
double precision, dimension(1:2,0:points) :: tmpt
double precision :: tmp
integer :: i,left,right,j
tmp1 = (/(y(i)*(h*dble(i))**2,i=0,points)/) !Première normalisation: intégrale 1d avec Jacobien radial 3d
call traptab(tmp1,h,points,tmp)
y = y/(tmp*4.d0*pi)

!normalisation du second moment: choisit dans la ligne de points fixes u**3*f(u*x) le paramètre u tel que le second moment soit normalisé à 1 (le 4eme moment 1d à 1/(4pi)). Soit a = \int_0^\infty x^4 f(x), alors u = \sqrt(4 \pi a).
tmp1 = (/(y(i)*(h*dble(i))**4,i=0,points)/) !intégrale 1d avec Jacobien radial 3d; cherche u (ou lambda)
call traptab(tmp1,h,points,tmp)
if (tmp<0) write(*,*) "ATTENTION: second moment < 0 (queue négative)."
tmp = dsqrt(4.d0*pi*dabs(tmp)) !tmp = valeur de u (ou de lambda)
l = tmp
if (tmp < 1.d0) then !Etirement vers la droite. Idée: construit le tableau f(u*x) puis interpole lin. les points
	tmpt(1,0:points) = (/(h*dble(i)/tmp,i=0,points)/) !valeurs de x étiré
	tmpt(2,0:points) = y !valeurs de f(x) correspondant
	tmp1(0) = y(0) !l'étirement n'affecte pas l'origine
	j=0
	do i=1,points !trouver h*i dans I=[h*j/tmp; h*(j+1)/tmp] et interpolation lin. pour tmp1(i) dans I
		recherche: do
			if ((dble(j) <= tmp*dble(i)) .and. (dble(j+1) >= tmp*dble(i))) then
				tmp1(i) = IL(tmpt(1,j),tmpt(1,j+1),tmpt(2,j),tmpt(2,j+1),dble(i)*h)
				exit recherche
			end if
			j=j+1
		end do recherche
	end do
	y = tmp1*tmp**3
end if
if (tmp > 1.d0) then !compression vers la gauche
	do i=1,nint(dble(points)/tmp)-1 !regarde les pts doubles (=envoyés s/même nouveau pt) pour en faire la moyenne
		left = nint(dble(i)*tmp) - nint(dble(i-1)*tmp) !nombre de tels points -1 vers la gauche
		right = nint(dble(i+1)*tmp) - nint(dble(i)*tmp) !nombre de tels points -1 vers la droite
		tmp1(i) = sum( y( nint(dble(i)*tmp)-left+1:nint(dble(i)*tmp)+right-1 ) )/dble(left+right-1) !moyenne
	end do
	right = nint(dble(1)*tmp) !cas i=0: seulement right
	tmp1(0) = sum(y(0:right-1))/dble(right)
	left = points - nint(dble(nint(dble(points)/tmp)-1)*tmp) !cas i=nint(dble(points)/tmp): seulement left
	tmp1(nint(dble(points)/tmp)) = sum(y(points-left+1:points))/dble(left)
	!points "manquants" = dernier point calculé (queue constante). A AMELIORER: cf notes.
	if (nint(dble(points)/tmp) < points) tmp1(nint(dble(points)/tmp)+1:points) = tmp1(nint(dble(points)/tmp))
	y = tmp1*tmp**3
end if

tmp1 = (/(y(i)*(h*dble(i))**2,i=0,points)/) !Seconde normalisation: intégrale 1d avec Jacobien radial 3d
call traptab(tmp1,h,points,tmp)
y = y/(tmp*4.d0*pi)
return
end subroutine normalise

!------------------------------------------------------------------------
!Ecriture de la solution dans le fichier "ann-eq.txt"
!  y       : tableau de dimension (n+1) contenant les valeurs y(xi)
!  h       : largeur du pas d'intégration
subroutine savedata(y,h)
implicit none
double precision, intent(in), dimension(0:points) :: y
double precision, intent(in) :: h
integer :: i
open(unit=1, file='ann-eq.txt')
do i=0,points
	write(1,*) h*dble(i),y(i)
end do
close(1)
return
end subroutine savedata

!------------------------------------------------------------------------
!Calcul des moments
!  y       : tableau de dimension (n+1) contenant les valeurs y(xi)
!  h       : largeur du pas d'intégration
!  M12     : moment <c_{12} c^2>
!  M10     : moment <c_{12}>
!  M02     : moment <c^2>
!  alpha   : rapport M12/(M10*M02)
subroutine moments(y,h,M12,M10,M02,alpha)
implicit none
double precision, intent(in), dimension(0:points) :: y
double precision, intent(in) :: h
double precision, intent(inout) :: M12,M10,M02,alpha
double precision, dimension(0:points) :: tmp1,tmp2
double precision :: tmp
integer :: i,j
!calcul de M12
do i=0,points
	tmp1=(/(y(i)*y(j)*(h*dble(i))**3*(h*dble(j))*((h*dble(i)+h*dble(j))**3-dabs(h*dble(i)-h*dble(j))**3),j=0,points)/)
	call traptab(tmp1,h,points,tmp2(i)) !seconde intégrale double [0,infty[: Gi = \int dj Fij
end do
call traptab(tmp2,h,points,tmp) !première intégrale double [0,infty[: H = \int di Gi
M12 = 8.d0*pi**2*tmp/3.d0

!calcul de M10
do i=0,points
	tmp1=(/(y(i)*y(j)*(h*dble(i))*(h*dble(j))*((h*dble(i)+h*dble(j))**3-dabs(h*dble(i)-h*dble(j))**3),j=0,points)/)
	call traptab(tmp1,h,points,tmp2(i)) !seconde intégrale double [0,infty[: Gi = \int dj Fij
end do
call traptab(tmp2,h,points,tmp) !première intégrale double [0,infty[: H = \int di Gi
M10 = 8.d0*pi**2*tmp/3.d0

!calcul de M02
tmp1 = (/(y(i)*(h*dble(i))**4,i=0,points)/) !construit f(c)*c^4
call traptab(tmp1,h,points,tmp) !première intégrale [0,infty[
M02 = tmp
tmp1 = (/(y(i)*(h*dble(i))**2,i=0,points)/) !construit f(c)*c^2
call traptab(tmp1,h,points,tmp) !seconde intégrale [0,infty[
M02 = M02*tmp*16.d0*pi**2

!calcul de alpha
alpha = M12/(M10*M02)
return
end subroutine moments

end module fonctions


!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------


program anneq
use fonctions
implicit none
double precision, dimension(0:points) :: f2,f1,f2b,tmp1,tmp2,tmp3
double precision :: h,alpha,M12,M10,M02,err,m2
integer :: i,j,k,it,largeur,one

h = (infty-0.d0)/dble(points) !doit être le même pour toutes les intégrations pour calculer les mêmes points
one = int(1.d0/h) !point le plus proche de 1 tout en restant <= 1.
largeur = 2*(one - nint(centrage/h)) !largeur totale du recollement en points
do i=0,points !condition initiale
	f1(i) = (3.d0/(2.d0*pi))**(3.d0/2.d0)*dexp(-3.d0/2.d0*(h*dble(i))**2)
end do
call normalise(f1,h,m2)
call savedata(f1,h)

do it=1,itmax !itérations
	call moments(f1,h,M12,M10,M02,alpha)

	!calcul du premier terme, petits c
	do i=0,one !valeurs de c
		tmp1(0:one) = (/(f1(nint(dble(i)*dble(j)*h))*(h*dble(j))**2,j=0,one)/) !choix arbitraire: point le + proche
!Interpolation linéaire: semble moins stable avec peu de points, Même résultats avec beaucoup de points.
!		tmp1(0) = 0.d0
!		tmp1(1:one) = (/( IL(dble(int(dble(i*j)*h**2)),dble(int(dble(i*j)*h**2))+1.d0,f1(int(dble(i*j)*h)),&
!												&f1(int(dble(i*j)*h)+1),dble(i*j)*h**2)*(h*dble(j))**2,j=1,one )/)
		call traptab(tmp1,h,one,tmp2(i))
	end do
	f2 = -tmp2
	
	!calcul du premier terme, grands c: intégration de [c,infty[
	tmp1 = (/(f1(i)*(h*dble(i))**2,i=0,points)/) !construit f(c)*c^2
	do i=one-largeur,points
		tmp2(0:points-i) = tmp1(i:points)
		call traptab(tmp2,h,points-i,tmp3(i))
	end do
	f2b = tmp3
	
	!calcul du second terme, petits c
	do i=0,one !valeurs de c
		do j=0,one
			tmp1 = (/(f1(nint(dble(i)*dble(j)*h))*f1(nint(dble(i)*dble(k)*h))*(h*dble(j))*(h*dble(k))*&
							&( (h*dble(j)+h*dble(k))**3 - dabs(h*dble(j)-h*dble(k))**3 ),k=0,points)/) !choix arbitraire: nint
!Interpolation linéaire: semble moins stable avec peu de points, Même résultats avec beaucoup de points.
!			tmp1 = (/( IL(dble(int(dble(i*j)*h**2)),dble(int(dble(i*j)*h**2))+1.d0,f1(int(dble(i*j)*h)),&
!										&f1(int(dble(i*j)*h)+1),dble(i*j)*h**2)&
!						   &*IL(dble(int(dble(i*k)*h**2)),dble(int(dble(i*k)*h**2))+1.d0,f1(int(dble(i*k)*h)),&
!										&f1(int(dble(i*k)*h)+1),dble(i*k)*h**2)&
!							 &*(h*dble(j))*(h*dble(k))*( (h*dble(j)+h*dble(k))**3 - dabs(h*dble(j)-h*dble(k))**3 ),k=0,points)/)
			call traptab(tmp1,h,points,tmp2(j)) !intégrale double [0,infty[: Gj = \int dk Fjk
		end do
		call traptab(tmp2,h,one,tmp3(i)) !intégration de [0,1] sur x: H(c_i) = \int dj Gj
		tmp3(i) = tmp3(i)*(h*dble(i))**2
	end do
	f2 = f2 + tmp3*2.d0*pi/(3.d0*M10)

	!calcul du second terme, grands c
	do i=one-largeur,points
		tmp1=(/(f1(i)*f1(j)*(h*dble(i))*(h*dble(j))*((h*dble(i)+h*dble(j))**3-dabs(h*dble(i)-h*dble(j))**3),j=0,points)/)
		call traptab(tmp1,h,points,tmp2(i)) !seconde intégrale double [0,infty[: Gi = \int dj Fij
	end do
	do i=one-largeur,points
		tmp1(0:points-i) = tmp2(i:points) 
		call traptab(tmp1,h,points-i,tmp3(i)) !première intégrale double [x,infty[: H = \int_x^infty di Gi
	end do
	f2b = f2b - tmp3*2.d0*pi/(3.d0*M10)

	!calcul de f_{n+1}
	f2 = 2.d0/(1.d0 - alpha)*f2
	f2b(one-largeur:points) = (/(2.d0/(1.d0 - alpha)*f2b(i)/(h*dble(i))**3,i=one-largeur,points)/)

	!recollement
	f2(one-largeur:one) = (f2(one-largeur:one) + f2b(one-largeur:one))/2.d0
	f2(one+1:points) = f2b(one+1:points)

	!calcul des erreurs, affichage, et de f_{n+1} avec lambda
	call normalise(f2,h,m2)
	err = sum(dabs(f1-f2))
	f1 = (1.d0 - lambda)*f1 + lambda*f2
	call normalise(f1,h,m2)
	call savedata(f1,h)
	write(*,'(" Iteration: ",i4,", err.abs.tot.: ",e8.2,", err.abs./point: ",e8.2,", l[1]: ",f5.2)') it,err,err/dble(points),m2 
end do

write(*,*) "Fin du programme."
read(*,*)
end program anneq